# ------------------------------------------------------------------------------#
# Description: Create:
#   Table A.4: Peer effects - rain garden eligible
#   Table A.5: Peer effects heterogeneity by type of GSI - rain garden eligible
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(data.table)
library(fixest)
library(modelsummary)

options(scipen = 999)

# Load and merge base data -----------------------------------------------------
load("data/adoptions/adoptions.RData")
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale    := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale     := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale     := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale   := peer.adopt.both.100 / 100]

dat <- dat[year > 2010]
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Table A.4: Peer effects - rain garden eligible ------------------------------------------------

dat[, peer.adopt.any.100.ols := peer.adopt.any.100]

m1 <- feols(adopt.any * 100 ~ peers.100 + peer.adopt.any.100.ols | year + blkgrp,
            cluster = 'blkgrp',
            data = dat[elig.rg == 1 & post.rainwise == 0])

m2 <- feols(peer.adopt.any.100 ~ peers.100 + elig.peers.cs.100 + elig.peers.rg.100 | year + blkgrp,
            cluster = 'blkgrp',
            data = dat[elig.rg == 1 & post.rainwise == 0])

m3 <- feols(adopt.any * 100 ~ peers.100 | year + blkgrp |
              peer.adopt.any.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
            cluster = 'blkgrp',
            data = dat[elig.rg == 1 & post.rainwise == 0])

m4 <- feols(adopt.any * 100 ~ peers.100 | year + blkgrp |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = 'blkgrp',
            data = dat[elig.rg == 1 & post.rainwise == 0])

varnames <- c(
  'year' = 'Year',
  'blkgrp' = 'Block Group',
  'peer.adopt.any.100.ols' = 'Peer Adoptions',
  'elig.peers.cs.100' = '\\# Eligible Peers (CS)',
  'elig.peers.rg.100' = '\\# Eligible Peers (RG)',
  'peer.adopt.any.100' = '$\\widehat{\\text{Peer Adoptions}}$',
  'peer.adopt.any.100.scale' = '$\\widehat{\\text{Peer Adoptions}}$'
)

# Load diagnostics
load("data/iv/iv_diag_table3_4_rg.RData")

eff.f.mult  <- round(iv_diag_table3_4_rg[['table3.iv.mult']]$F_stat[4], 2)
eff.f.sing  <- round(iv_diag_table3_4_rg[['table3.iv']]$F_stat[4], 2)
ar.ci.mult  <- paste0("[", round(iv_diag_table3_4_rg[['table3.iv.mult']]$AR$ci[1], 2), ", ",
                      round(iv_diag_table3_4_rg[['table3.iv']]$AR$ci[2], 2), "]")
ar.ci.sing  <- paste0("[", round(iv_diag_table3_4_rg[['table3.iv']]$AR$ci[1], 2), ", ",
                      round(iv_diag_table3_4_rg[['table3.iv']]$AR$ci[2], 2), "]")
tf.ci.sing  <- paste0("[", round(iv_diag_table3_4_rg[['table3.iv']]$tF[6], 2), ", ",
                      round(iv_diag_table3_4_rg[['table3.iv']]$tF[7], 2), "]")
tf.cf.sing  <- round(iv_diag_table3_4_rg[['table3.iv']]$tF[2], 2)

# Preview Table A.4
etable(m1, m2, m3, m4,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       drop = c("peers.100"),
       fitstat = c("n", "sargan", "wh"),
       extralines = list(
         '_F' = c("", "", eff.f.mult, eff.f.sing),
         '_AR CI' = c("", "", ar.ci.mult, ar.ci.sing),
         '_tF cF' = c("", "", "", tf.cf.sing),
         '_tF CI' = c("", "", "", tf.ci.sing)
       ))

# Save Table A.4
style_tex <- style.tex(
  tablefoot.value = "",
  model.title = "",
  depvar.title = "",
  var.title = "\\midrule",
  stats.title = "\\midrule"
)

etable(m1, m2, m3, m4, tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       drop = c("peers.100"),
       fitstat = c("n", "sargan", "wh"),
       headers = list("_ " = c("OLS", "First Stage", "IV", "IV")),
       extralines = list(
         '_F' = c("", "", eff.f.mult, eff.f.sing),
         '_AR CI' = c("", "", ar.ci.mult, ar.ci.sing),
         '_tF cF' = c("", "", "", tf.cf.sing),
         '_tF CI' = c("", "", "", tf.ci.sing)
       ),
       file = "output/tables/table_a4.tex", replace = TRUE)

rm(m1, m2, m3, m4, ar.ci.mult, ar.ci.sing, eff.f.mult, eff.f.sing, tf.cf.sing, tf.ci.sing)

# Table A.5: Peer effects heterogeneity by type of GSI - RG eligible ------------------------------------------------

base <- feols(adopt.any ~ peers.100 | year + blkgrp |
                peer.adopt.any.100.scale ~ elig.peers.avg.100,
              cluster = "blkgrp",
              data = dat[elig.rg == 1 & post.rainwise == 0])

m1 <- feols(adopt.any.rg ~ peers.100 | year + blkgrp |
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.rg == 1 & post.rainwise == 0])

m2 <- feols(adopt.cs ~ peers.100 | year + blkgrp |
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.rg == 1 & post.rainwise == 0])

m3 <- feols(adopt.rg ~ peers.100 | year + blkgrp |
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.rg == 1 & post.rainwise == 0])

m4 <- feols(adopt.both ~ peers.100 | year + blkgrp |
              peer.adopt.any.100.scale ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.rg == 1 & post.rainwise == 0])

eff.f.1 <- round(iv_diag_table3_4_rg[['table3.iv']]$F_stat[4], 2)
eff.f.2 <- round(iv_diag_table3_4_rg[['table4.any.rg']]$F_stat[4], 2)
eff.f.3 <- round(iv_diag_table3_4_rg[['table4.cs']]$F_stat[4], 2)
eff.f.4 <- round(iv_diag_table3_4_rg[['table4.rg']]$F_stat[4], 2)
eff.f.5 <- round(iv_diag_table3_4_rg[['table4.both']]$F_stat[4], 2)

ar.ci.1 <- paste0("[", round(iv_diag_table3_4_rg[['table3.iv']]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4_rg[['table3.iv']]$AR$ci[2], 2), "]")
ar.ci.2 <- paste0("[", round(iv_diag_table3_4_rg[['table4.any.rg']]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4_rg[['table4.any.rg']]$AR$ci[2], 2), "]")
ar.ci.3 <- paste0("[", round(iv_diag_table3_4_rg[['table4.cs']]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4_rg[['table4.cs']]$AR$ci[2], 2), "]")
ar.ci.4 <- paste0("[", round(iv_diag_table3_4_rg[['table4.rg']]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4_rg[['table4.rg']]$AR$ci[2], 2), "]")
ar.ci.5 <- paste0("[", round(iv_diag_table3_4_rg[['table4.both']]$AR$ci[1], 2), ", ",
                  round(iv_diag_table3_4_rg[['table4.both']]$AR$ci[2], 2), "]")

# Preview Table A.5
etable(base, m1, m2, m3, m4,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       headers = list("_ " = c("Any", "Rain Garden", "Only Cistern", "Only Rain Garden", "Both")),
       fitstat = c("n"),
       extralines = list(
         '_F' = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5),
         '_AR CI' = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5)
       ))

# Save Table A.5
etable(base, m1, m2, m3, m4, tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       headers = list("_ " = c("Any", "Rain Garden", "Only Cistern", "Only Rain Garden", "Both")),
       fitstat = c("n"),
       extralines = list(
         '_F' = c(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5),
         '_AR CI' = c(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5)
       ),
       file = "output/tables/table_a5.tex", replace = TRUE)
